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For a two-dimensional system of coupled oscillators, 
the spectra of reduced density matrices can be obtained 
analytically. This provides an example where the fea- 
tures of these quantities, which are of central impor- 
tance in numerical studies using the DMRG method, 
can be seen. 

The density-matrix renormalization group 
method (DMRG) [Q has brought an enormous 
progress for the study of one-dimensional quantum 
systems and related classical problems Conse- 
quently, this numerical approach has also been ;ip- 
plied to two-dimensional quantum problems 
The situation there, however, appears to be much 
less favourable. 

The reason for this has to be sought in the prop- 
erties of the density matrices which are used to se- 
lect an optimal reduced basis in the Hilbert space. 
The essential quantity is the distribution of their 
eigenvalues. In one dimension, one usually finds 
a rapid decay, so that a relatively small number 
of states is sufficient to give very good results. 
This basically exponential decay can be derived 
explicitely for non-critical integrable models [^|J|]. 
In non-integrable cases, the spectra are less regular 
but still have similar features [0J|]. 

The situation in two dimensions has been dis- 
cussed in some detail for free fermions |^ and for 
the transverse Ising model It was found that, 
if one couples one-dimensional chains to form lad- 
ders, the number m of states one needs to maintain 
a certain accuracy, grows exponentially with the 
width M of the system. This was derived either 
from the limit of non-interacting chains, or from 
numerical calculations. The spectra themselves, 
however, have not been discussed so far, although 
they are at the core of the problem. It therefore 
seems worth while to treat an example, where one 
can give explicit results. 

This is possible for a system of coupled harmonic 
oscillators, which is integrable in any number of di- 
mensions. This problem was studied recently for 
the case of a linear chain Q and it was shown that 
the ground-state density matrices, either for one 



site or for half of the system, are exponentials of 
bosonic operators. This is a consequence of the 
Gaussian form of the ground state and holds quite 
generally. The problem is only to determine the 
bosonic eigenvalues. This can be done either nu- 
merically for a small system, or analytically in the 
thermodynamic limit. 

To be specific, consider the system described by 
the Hamiltonian 



(1) 



where Ui is the coordinate of the i-th oscillator 
and Wo its frequency. The masses are all equal 
to unity and the oscillators are coupled by springs 
of strenght fcy . Transforming H to normal coordi- 
nates, one can write down the ground state imme- 
diately. In terms of the original coordinates it has 
the form 



exp(-i^ A^UiUj) 



(2) 



The total density matrix is then 10 >< 0|. By 
integrating out part of the coordinates, one obtains 
reduced density matrices which have the diagonal 
form 



C exp (- 
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(3) 



with bosonic operators The eigenvalues 

Ej follow from a matrix which is obtained from 
Aij and has a dimension equal to the number 
of kept sites. One divides Aij into four subma- 
trices a^^ , a^^ , a^^ , a^^ , according to whether the 
sites i and j are kept or not. Then the matrix 
a^^[a^'^{a'^'^)~^a'^^]~^ has eigenvalues cosh2(ej72). 
This can be shown by a straightforward general- 
ization of the approach in section 2 of Q. In this 
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way, density-matrix spectra can be calculated nu- 
merically for an arbitrary assembly of coupled os- 
cillators. 

For a large system, however, the situation sim- 
plifies. In 1^] it was shown that, for a chain with 
nearest-neighbour coupling k and oscillator fre- 
quency ujQ = I — k, the Ej for half of the system 
are, in the thermodynamic limit, given by 

e, = (2j-l)e, J -1,2,3,... (4) 

where 

e = nl{k')/l{k) (5) 

Here I{k) is the complete elliptic integral of the 
first kind and k' = \/\ — k"^ . The result is also valid 
for finite systems if the size is large compared with 
the correlation length. The Sj for smaller systems 
are still similar, but there are deviations from 
which increase for larger values of j . 

Now consider a two-dimensional square lattice 
of oscillators with nearest-neighbour couplings 
and ky in the two directions. This can be reduced 
to a one-dimensional problem by first introducing 
normal coordinates in the columns. The corre- 
sponding normal frequencies are 

ujlyqf = bjl^2ky{\- co^q) (6) 

where the vertical momenta for open boundary 
conditions at the ends and M sites are given by 

fry 

q = —7T, m = 0,l,2,...,(Af-l) (7) 

If one now couples the columns, the different 
momenta do not mix, and for each value of g a 
horizontal chain of the form P) results where the 
oscillator frequency is now i^(q) and the coupling 
kx- For the density matrix of the half-system, this 
leads to the spectrum (0),(||) with the parameter 
k = k{q) determined from the relation kx/uj{q) — 
fc/(l — k) or explicitely 

k = k^/{k^+Lo{q)) (8) 

In this way, an analytic expression for the spec- 
trum is obtained. For each j, one has a band of 
M eigenvalues ej{q) due to the transverse exten- 
sion of the system. This reflects the corresponding 
interface between the two parts into which the sys- 
tem is divided. The dispersion of the vibrational 



modes in the vertical direction also determines the 
dispersion of the e-band via (H). In particular, a 
large Lo{q) also leads to a large value of e. 
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FIG. 1. Bosonic single-particle eigenvalues ej{q), 
c.f. Equ.Q, for half of a 10 x 10 system, arranged in 
ascending order, for uio = kx = 1.0 and four values of 
the coupling ky. 



Such spectra, calculated numerically for a 10 x 10 
lattice, are show in Fig.|l|. Plotted are the ej{q) 
arranged according to their magnitude for differ- 
ent values of the transverse coupling ky . For non- 
interacting horizontal chains, one has a sequence 
of plateaus. Turning on ky{> 0), the eigenvalues 
increase except for q = and form real bands. At 
the lower end, the stair-like structure still persists, 
while for larger values of j the bands are spread 
more due to the factor {2j — 1) and overlap even- 
tually. After the proper ordering of the e, a con- 
tinuous curve emerges. It corresponds roughly to 
a linear relation of the form 

£„ = An (9) 

with integer n and A = 2e{q = 0)/M inversely 
proportional to the width M. 

The actual eigenvalues Wn of p are obtained by 
specifying the occupation numbers of the bosonic 
single-particle levels £j{q). This leads to an in- 
creasingly larger number of possibilities as more 
e are involved, i.e. for smaller Wn- The final re- 
sult is shown in Fig.^ for the same parameters as 
in Fig.0. One can see that the stair-like structure 
persists also in the Wn for small ky, although the 
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plateaus are much longer and given by combina- 
torical factors. For larger ky, rather smooth curves 
arise which drop increasingly faster. In all cases, 
there is a rapid initial decay followed by a slower 
decrease for larger n. Following one can derive 
an asymptotic formula from (O) which reads 



exp{-(A/(2 7rV3))ln2n} 



(10) 



and which is obeyed reasonably well by the curves. 
Due to the slow decay, the truncation error when 
cutting off the spectrum also decreases slowly. Af- 
ter n = 100, 500 and 1000 it is approximately 10~^, 
10"''^ and 10^*, respectively, if kx = ky = loq = 1.0. 



k =0.0 

k'=0.1 

k'=1.0 

k=10.0 



such that the Hamiltonian expressed in Fermi op- 
erators has the same structure as (|^) expressed in 
Bose operators. In this case, the combinatorical 
possibilities are reduced, but this leads only to a 
change A ^ 2A in (p^. One should also note that 
we have treated a non-critical system where the 
situation is in general more favourable. One could 
extend the considerations to three dimensions, in 
which case one has two momenta for the transverse 
directions and therefore an even larger number of 
£-values in each band. 
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FIG. 3. Density-matrix eigenvalues Wn for systems 



of different width M and loq = kx — ky 



1.0 



FIG. 2. Density-matrix eigenvalues w„, arranged 
in decreasing order, obtained from the £j{q) in Fig.|l| 
and for the same parameters 



The dependence of the ui^-spectrum on the 
width M is shown in Fig.^for the case k^ = ky = \. 
One can see how the curves drop more and more 
slowly as M increases, in accordance with and 
( |l0| ) and the decrease of A with M . 

These results confirm that the situation worsens 
as the system becomes more two-dimensional. The 
faster initial decay before the onset of the combina- 
torical effects helps in numerical calculations. Also 
the interaction helps here to some extent since the 
e- values increase with ky , but this does not remove 
the basic dependence in the exponent. The 
same features are found if one assumes that one is 
dealing with fermionic operators in . This would 
correspond to a fermionic system with pair terms 



Returning briefly to one dimension, we would 
like to mention that our model also shows the ori- 
gin of the different DMRG performance for chains 
and rings ||l| very clearly. This problem is, in fact, 
closely related to those discussed above. If one cal- 
culates the density-matrix spectrum for a half-ring, 
one finds the same -values at the lower end as for 
the half-chain, but each value appears twice. When 
plotted, this leads to a structure as in Figj^ with 
steps of length two. The reason lies in the form of 
the eigenstates of p which, for small Sj, are con- 
centrated near the boundary between system and 
environment. This feature, which was conjectured 
before 11|, can be seen explicitely here and is 
illustrated in Fig.^. The effect is known from the 
closely related corner transfer matrix of the mas- 
sive Gaussian model |l^. For a half-ring, which 
has two points of contact, one then finds two such 
sets of states which are approximately independent 
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of each other for smaU Ej. Therefore, p ^ Pl ■ Pb. 
where pL and pa are density matrices for only a 
left or only a right bondary. Thus the situation 
is the same as for a ladder consisting of two only 
weakly interacting chains. 




FIG. 4. Density-matrix eigenstates for the left part 
of a chain of 32 sites. Shown are the amplitudes as a 
function of the position for the lowest three Sj-values 
for LOO = k = 0.5. 
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Coming back to the two-dimensional problem, 
the spectra found here explain, in a very direct 
way, the difficulties of the DMRG in this case. To 
apply the method, one should use as many symme- 
tries as possible [Q. However, to treat really large 
systems, a procedure which avoids the extended 
interfaces between the parts of the system would 
be necessary. Whether the momentum-space ap- 
proach of Ref. U] can help here, is not yet clear. 
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